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Abstract.  The  transport  coefficient  coupling  heat  flux  and  density  gradient  in  a  granular  gas  is  measured  by  taking  advantage 
of  the  existence  of  a  minimum  in  the  temperature  profile  of  an  open  vibrated  granular  medium.  This  temperature  inversion  is 
closely  related  to  the  existence  of  the  new  transport  coefficient.  Particle  simulations  using  the  Direct  Simulation  Monte  Carlo 
method  will  be  used  to  compute  the  transport  coefficient,  and  the  results  will  be  compared  with  theoretical  predictions  derived 
from  the  Boltzmann  equation.  Finally,  the  accuracy  of  a  boundary  condition  requiring  the  hydrodynamic  heat  flux  to  vanish 
for  infinite  heights  will  be  discussed. 


INTRODUCTION 

Steady  states  of  granular  materials  are  obtained  when  energy  is  supplied  to  the  system  in  order  to  compensate  the 
energy  loss  in  collisions.  One  of  the  easiest  ways  to  add  energy  to  a  system  in  experiments  is  to  do  so  through  a 
vibrating  wall.  If  the  vibration  is  strong  enough,  the  system  will  be  fluidized,  and  one  can  expect  its  behavior  to 
be  described  by  the  extension  of  the  Navier-Stokes  equations  to  granular  systems.  In  this  work,  we  will  use  the 
hydrodynamical  description  to  study  the  steady  state  of  an  open,  vibrated,  granular  system  in  presence  of  gravity. 
The  granular  fluid  will  be  modelled  as  a  system  of  inelastic  hard  particles,  whose  collisions  are  characterized  by  a 
constant  coefficient  of  normal  restitution  a.  Besides,  we  will  be  interested  in  the  dilute  limit,  when  the  Boltzmann 
equation  applies. 

A  distinctive  feature  of  granular  gases  as  compared  to  molecular  (elastic)  ones,  is  that  the  expression  of  the  heat 
flux  has  to  be  generalized  by  including  a  new  term  coupling  the  heat  flux  and  the  density  gradient.  This  implies  the 
introduction  of  a  new  coefficient,  the  diffusive  heat  conductivity  /i,  that  vanishes  in  the  elastic  limit.  This  coupling  has 
been  derived  by  kinetic  theory  methods  [1,  2,  3],  and  its  consequences  confirmed  in  computer  simulations  [4,  5,  6]. 
For  the  particular  case  of  a  vibrated  granular  gas  in  the  presence  of  gravity,  this  term  implies  a  peculiar  behavior  of 
the  temperature  profile,  that  increases  with  height  after  a  minimum.  The  existence  of  the  minimum  was  first  derived 
from  the  hydrodynamic  equations  [5],  and  was  confirmed  by  computer  simulations  [5,  7,  8],  and  also  in  experiments 
[9].  Here,  we  will  show  that  the  value  of  the  diffusive  heat  conductivity  /i  can  be  obtained  from  the  behavior  of  the 
system  at  the  temperature  minimum.  Then,  the  value  of  the  transport  coefficient  will  be  compared  with  the  theoretical 
prediction  from  the  Boltzmann  equation  derived  in  [2],  Finally,  the  boundary  conditions  to  be  used  when  solving  the 
hydrodynamic  equations  will  be  also  discussed. 


HYDRODYNAMIC  DESCRIPTION 

Let  us  consider  a  system  of  N  inelastic  hard  spheres  ( d  =  3)  or  disks  (d  =  2)  of  mass  m  and  diameter  cr,  in  presence  of  a 
gravitational  field  g  =  —  gez,  where  g  is  a  positive  constant  and  ez  a  unit  vector  in  the  Z  direction.  In  the  hydrodynamic 
description,  the  state  of  the  system  is  completely  specified  by  the  local  number  of  particles  density,  n(r,t),  the  velocity 
flow  u(r,f),  and  the  temperature,  T(r.t).  For  a  dilute  gas,  the  evolution  of  these  quantities  is  given  by  the  extension  to 
the  inelastic  case  of  the  Navier-Stokes  equations  [2,  10], 

|^  +  V-(m i)  =  0,  (1) 
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+  u  •  Vu  +  —V  •  ^  —  g  =  0, 
dt  nm 

d  T  2 

—  +  u •  VT  +  — [<?  :  Vu  +  V •  q]  +  =  0 . 

dt  ankB 

Here,  k/>,  is  the  Boltzmann  constant,  is  the  pressure  tensor. 


1  =  pj?  -  p 


(Vu)  +  (Vu)+--J*rV-u 


(3) 

(4) 


where  p  =  nkuT  is  the  hydrostatic  pressure,  ■  '/  the  identity  tensor,  and  rj  the  shear  viscosity  coefficient,  q  appearing 
in  Eq.  (3)  is  the  heat  flux, 

q  =  —kVT  —  pVn ,  (5) 

with  K  the  heat  conductivity  coefficient,  and  p  the  diffusive  heat  conductivity.  Finally,  £  is  the  cooling  rate  associated 
to  the  energy  dissipation  in  collisions.  The  expression  for  these  coefficients  reads: 


r]  =  Tj*(a)rjo(T),  K  =  K*(a)Ko(T),  p  =  p*(a)p0  (T),  (6) 

c^c(0)  =  n«)f,  (?) 

with  rjo,  K"o  the  Boltzmann  elastic  values  of  the  viscosity  and  heat  conductivity,  po  =  TKo/n,  while  rj*(a),  K*(a ), 
p*(oc),  and  C*(°0  are  dimensionless  functions  of  the  coefficient  of  restitution.  Although  their  expressions  will  not  be 
given  here  (they  can  be  found  in  [2,  10]),  it  is  important  to  remember  that  in  the  elastic  limit  a  — >  1,  rj*  and  K*  tend  to 
unity,  while  p*  and  vanish.  Besides,  both  rjo  and  K'o  are  proportional  to  7’ 1  2. 

Let  us  consider  that  energy  is  supplied  to  the  system  through  a  vibrating  wall  of  size  S  located  at  z  =  0.  Although 
the  details  of  the  wall  movement  are  not  very  important  here,  we  will  consider  that  the  vibration  amplitude  is  small 
enough  as  to  approximate  the  position  of  the  wall  as  fixed.  Also,  the  wall  moves  with  a  sawtooth  profile,  so  particles 
find  it  moving  upwards  with  a  characteristic  velocity.  When  the  energy  supplied  by  the  wall  compensates  the  one  lost 
in  collisions,  the  system  reaches  a  steady  state.  Besides,  and  because  of  the  symmetry  of  the  problem,  we  can  expect 
that,  once  in  this  state,  there  will  be  gradients  only  in  the  Z  direction.  When  Eqs.  (l)-(3)  are  particularized  for  this 
state,  they  take  the  form 


dp 

j-  =  —rung, 

(8) 

*  d(!Z  +  TC(0)  =  0. 

ankB  dz 

(9) 

Besides,  Eq.  (8)  implies  that  the  heat  flux  in  this  case  is 

qz(z)  =  -(K*-p*)K0^  +  p*Ko'^-. 

dz  kg 

(10) 

In  order  to  solve  these  equations,  it  is  convenient  to  introduce  the  new  dimensionless  length  scale  <5  as: 


where  A(z)  is  the  local  mean  free  path  for  hard  disks  or  spheres,  A(z)  =  [Cncr^-1]-1,  with  C  =  2\/2  for  d  =  2  and 
C  =  7T%/2  for  d  =  3.  Let  us  notice  that,  when  z  — >  ^  — *  0,  while  E,  takes  its  maximum  value,  E,q  =  y/aCo‘l^1Nz,  with 

Nz  =  N/S,  at  z  =  0.  The  coefficient  a  appearing  in  Eq.  (11)  is  a  function  of  the  coefficient  of  restitution  and  reads 


32(d  —  \)nd~l  C  *(«) 

C2(c/  +  2)3r((7/2)2  K*(a)-p*(a)  ’ 


(12) 


vanishing,  therefore,  in  the  elastic  limit.  In  terms  of  the  new  scale,  the  solution  to  the  Navier-Stokes  equations  is  [5]: 


T1/2(Z)  =  rVlAIv(Z)+BKv(Z)\, 


(13) 
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FIGURE  1.  Temperature  and  density  profiles  for  a  system  of  hard  spheres  with  a  =  0.925,  =  1.92.  The  symbols  are  from 

the  simulations,  while  the  lines  are  the  solution  to  the  hydrodynamic  equations,  with  the  arbitrary  constants  determined  from  the 
temperature  minimum. 


*($)  = 


mgt, 


l+2v 


(14) 


CkBad^ y/a(a)  [AIV^)+BKV^)}2  ’ 

where  Iv  and  Kv  are  the  modified  Bessel  functions  of  first  and  second  kind,  and  A  and  B  are  constants  to  be  determined 
from  the  boundary  conditions.  The  parameter  v  is 


v(a) 


/**(«) 

4[ic*(a)  -  ju*(a)] 


>0, 


(15) 


In  the  limit  |  — >  0  (z  — >  °°)  Iv  — >  0,  while  Kv  — >  °°  [1 1],  Nevertheless,  this  does  not  imply  that  the  constant  B  has  to  be 
taken  identically  equal  to  zero,  as  there  is  nothing  unphysical  in  the  fact  that  T  diverges  as  far  as  the  density  decreases 
fast  enough  as  to  guarantee  that  the  local  kinetic  energy  density  goes  to  zero  in  that  limit.  Besides,  it  is  clear  that 
the  hydrodynamic  description  will  not  be  valid  for  very  large  heights,  as  the  density  will  have  decayed  to  very  small 
values,  and  the  local  Knudsen  number  will  be  very  small.  Then,  we  will  keep  the  B  constant  different  from  zero. 

The  presence  of  the  term  proportional  to  Kv  in  the  temperature  profile  implies  that  it  has  a  minimum  located  at 
4  =  4n  given  by 

AIv+1($m)-BKv+l($m)  =  0,  (16) 

and  the  temperature  Tm  at  the  minimum  is 

Tm1  =  4,7 V  Wvi^m)  +BKv($m)} .  (17) 


Then,  if  hydrodynamics  is  valid  in  the  vicinity  of  the  temperature  minimum,  the  constants  A  and  B  can  be  determined 
from  the  measured  temperature  profiles.  It  has  been  also  shown  [5]  that  the  density  has  a  maximum  at  t,  =  that  is 
approximately  given  by  the  solution  of  the  equation: 

/v(4,)-24,/v+1(4,)  =  o,  (18) 

which  can  be  solved  numerically  for  each  value  of  a.  Let  us  just  comment  that  4,  takes  values  of  the  order  of  unity. 
Of  course,  in  order  to  observe  the  density  maximum  in  an  experiment  the  number  of  particles  in  the  system  has  to  be 
large  enough  so  that  |o  >  4,-  If  this  condition  is  not  fulfilled,  the  density  will  decay  monotonically  with  height. 

In  order  to  check  the  above  hydrodynamic  description,  computer  simulations  for  two  and  three  dimensional  systems 
by  using  the  Direct  Simulation  Monte  Carlo  method  (DSMC)  [12]  have  been  performed.  This  method  is  particularly 
suited  for  this  problem  as  we  are  interested  in  the  low  density  limit,  and  as  it  allows  to  exploit  the  symmetry  of  the 
system  in  the  transversal  direction.  In  Fig.  1  the  steady  density  and  temperature  profiles  for  a  system  of  hard  spheres 
have  been  plotted.  The  density  has  been  scaled  with  the  initial,  homogeneous  density,  while  the  temperature  is  scaled 
with  some  arbitrary  value.  Space  is  measured  in  units  of  the  initial,  homogeneous,  mean  free  path.  The  values  of  the 
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FIGURE  2.  Reduced  transport  coefficient  fl*  as  a  function  of  a  for  a  system  of  hard  spheres.  The  solid  line  is  the  theoretical 
prediction  derived  in  [2],  and  the  circles  the  results  from  the  simulation  of  a  vibrated  system.  The  squares  are  also  simulation  results 
but  using  a  Green-Kubo  expression  for  fl. 


parameters  of  the  simulations  were  a  =  0.925,  c,o  =  1.92  .  The  symbols  are  the  results  of  the  simulations:  as  predicted, 
the  density  shows  a  maximum,  while  the  temperature  exhibits  a  minimum,  increasing  from  there  on.  The  continuous 
lines  are  the  theoretical  prediction,  with  A  and  B  determined  from  the  temperature  minimum.  The  agreement  is  very 
good,  confirming  the  validity  of  hydrodynamics  in  the  temperature  minimum  region.  The  same  qualitative  results  are 
achieved  for  other  values  of  the  parameters  both  in  the  two  and  three  dimensional  cases,  as  far  as  the  coefficient  of 
restitution  is  not  too  low.  Due  to  an  intrinsic  coupling  between  gradients  and  inelasticity  in  this  steady  state,  retaining 
only  up  to  the  Navier-Stokes  order  in  the  Chapman-Enskog  expansion  may  be  not  enough  for  small  values  of  a. 


THE  HEAT  FLUX 

In  the  steady  state  of  a  granular  material,  the  heat  flux  is  given  by  Eq.  (11).  At  the  temperature  minimum,  as  the 
derivative  of  T  vanishes,  we  have 

=  (19) 

kb 

Then,  the  computation  of  the  heat  flux  at  the  temperature  minimum  provides  a  direct  measurement  of  the  scaled 
diffusive  heat  conductivity  coefficient,  ji " .  It  must  be  pointed  out  that  Eq.  (19)  only  requires  the  validity  of  the 
hydrodynamic  description,  and  no  additional  condition  has  to  be  introduced.  In  Fig.  2  we  have  plotted  the  coefficient 
/i *  as  a  function  of  a  for  a  system  of  hard  spheres.  The  solid  line  is  the  theoretical  prediction  derived  in  [2]  by  using  the 
Chpaman-Enskog  expansion  from  the  Boltzmann  equation.  The  circles  are  the  results  of  the  simulation  using  Eq.  (19). 
The  error  bars  are  obtained  by  giving  different  values  to  the  parameter  go  and  the  velocity  of  the  vibrating  wall.  The 
agreement  between  theory  and  simulation  is  quite  good,  even  at  the  lowest  values  of  a  investigated.  Nevertheless, 
it  must  be  said  that,  for  those  values,  the  shape  of  the  profiles  begin  to  show  discrepancies  from  the  theoretical 
prediction  derived  here.  This  is  the  reason  of  the  large  error  bars  for  the  lowest  values  of  a,  and  for  not  having 
considered  smaller  values  of  this  parameter.  Finally,  we  have  also  included  in  the  figure  the  results  of  independent 
DSMC  simulations  where  fi*  was  computed  by  means  of  Green-Kubo  expressions  derived  in  Ref.  [13]  (squares).  The 
agreement  is  again  very  good.  It  could  be  argued  that  it  is  not  surprising  to  obtain  a  good  agreement  between  the  results 
of  the  DSMC  simulation  and  a  theoretical  description  derived  from  the  Boltzmann  equation.  Nevertheless,  it  must  be 
remembered  the  the  theoretical  derivation  requires  given  hypothesis  and  approximations  (validity  of  hydrodynamics, 
gradient  expansions,  Sonine  expansion)  that  are  not  assumed  at  all  by  the  simulation  method. 

It  is  important  to  stress  that  the  fact  that  the  heat  flux  does  not  vanish  at  the  temperature  minimum  is  a  direct  proof 
of  the  existence  of  the  coupling  between  heat  flux  and  density  gradient.  Besides,  Fig.  2  shows  that,  although  fi*  — >  0 
in  the  elastic  limit,  its  contribution  cannot  be  neglected  as  a  goes  beyond  the  quasi-elastic  limit. 
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FIGURE  3.  Position  of  the  temperature  minimum  t,m  and  of  the  density  maximum  £„in  the  three  dimensional  case.  The  lines  are 
the  theoretical  prediction  (solid  line  for  ,  dot-dashed  line  for  The  symbols  are  the  results  of  the  DSMC  simulations. 


Eq.  (19)  provides  the  heat  flux  at  the  temperature  minimum.  A  general  expression  for  qz  valid  at  any  position  can 
be  obtained  when  the  solution  given  by  Eq.  (14)  is  substituted  in  Eq.  (11),  and  it  reads 


qz=A(K*-ll*) 


2mgKb 
kBT '/2 


1 V— ,(4) 


:  1— v 


(20) 


Let  us  consider  now  the  limit  z  — >  °°  or,  equivalently,  %  — >  0.  Taking  into  account  the  asymptotic  behavior  of  the 
modified  Bessel  functions,  it  follows  that,  for  very  large  heights  is  given  by: 


qz  =  A(K*-n*) 


ImgKo 


kBTx/2 

If  we  require  that  the  heat  flux  vanishes  for  very  large  heights,  we  get  the  relation 


B  2 

a  ~  r(v)r(i-v)  ’ 


(21) 


(22) 


i.e.  the  ratio  of  the  two  constants  is  given  by  a  function  of  a  alone.  Taking  into  account  the  behavior  of  the  T  function, 
/I*  =  0  (i.e.  v  =  0)  implies  B  =  0,  and  the  temperature  would  be  constant  for  large  heights.  If  we  introduce  the  above 
relation  into  the  equation  that  determines  the  position  of  the  temperature  minimum,  Eq.  (16),  we  get  a  closed  equation 
for  %m\ 

A'+l  (^m)  —  r(y)r(l  —  y-j^v+1  (4z«)  =  0  •  (23) 

Then,  the  position  of  the  temperature  minimum  in  a  vibrated  system  in  the  scaled  space  variable  is  a  function  only  of 
the  coefficient  of  normal  restitution,  independent  of  the  other  relevant  parameters  of  the  system  (number  of  particles, 
gravity  and  vibration  velocity). 

In  Fig.  3  the  position  of  the  temperature  minimum  < in  the  three  dimensional  case  is  plotted  as  a  function  of  a. 
The  symbols  are  the  results  of  the  DSMC  simulation,  while  the  solid  line  is  the  solution  of  Eq.  (23).  The  agreement 
is  very  good,  supporting  the  use  of  the  above  mentioned  boundary  condition  to  solve  the  hydrodynamic  equations.  It 
must  be  noticed  that  the  validity  of  this  boundary  condition  is  not  at  all  clear,  as  it  is  imposed  in  the  very  large  heights 
region,  i.e.,  once  the  density  has  decayed  to  very  small  values  and  the  hydrodynamic  description  is  not  valid.  Of  course, 
the  heat  flux  must  vanish  in  that  region,  but  the  question  is  whether  this  can  be  translated  into  an  effective  boundary 
condition  for  the  hydrodynamic  equations.  For  instance,  a  vanishing  heat  flux  implies  a  constant  temperature  gradient 
whose  value  is  directly  related  to  v.  But  identifying  the  asymptotic  region  where  this  linear  behavior  is  achieved  is 
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very  difficult  in  practical  applications,  because  of  the  failure  of  hydrodynamics  to  describe  the  upper  region  of  the 
system  [6],  Here  we  have  shown  that  the  condition  of  vanishing  flux  at  large  heights  translates  into  a  condition  inside 
the  hydrodynamic  region,  so  it  can  be  clearly  tested. 

We  have  also  included  in  Fig.  3  the  position  of  the  density  maximum,  E,n.  The  dot-dashed  line  is  the  theoretical 
prediction  given  by  the  solution  of  Eq.  (18),  while  the  triangles  are  the  results  of  the  simulations.  The  agreement  is 
again  quite  good,  although  some  discrepancies  appear  for  the  lowest  values  of  a  studied.  The  reason  for  this  might 
be  that  Eq.  (18)  is  not  exact  [5],  and  is  in  fact  obtained  for  v<  1,  while  for  a  =  0.7,  v  ~  0.14.  Nevertheless,  it  must 
be  noticed  the  weak  dependence  on  a  of  as  compared  to  £,rn.  In  fact,  for  the  values  of  the  coefficient  of  restitution 
considered  in  the  simulations  <^„  ~  1 . 

In  conclusion,  hydrodynamics  provides  a  very  useful  tool  to  study  non-homogeneous  steady  states  of  granular 
systems.  Nevertheless,  this  hydrodynamic  description  has  distinctive  characteristics  that  cannot  be  guessed  from  the 
one  of  molecular  fluids.  In  particular,  a  new  transport  coefficient,  the  diffusive  heat  conductivity,  has  to  be  introduced. 
In  this  work  we  have  shown  that  this  new  transport  coefficient  has  relevant  consequences  in  the  hydrodynamic  profiles, 
so  it  cannot  be  neglected  in  the  description  of  granular  flows. 
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